Modeling Extreme Events

HTML(showcode)

The study domain and the model domain setup

#
fontsize_tmp = fontsize
fontsize=20
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes, inset_axes
from mpl_toolkits.axes_grid1.inset_locator import mark_inset
topocmap = plt.cm.gist_earth
topocmap.set_bad('w', alpha=0)
#First the greater area map
with nc(os.path.join(geodataf,'topo_2.nc')) as topof:
    Gtopo=np.ma.masked_equal(topof.variables['topo'][:],0)
    Gtopolon=topof.variables['longitude'][:]
    Gtopolat=topof.variables['latitude'][:]
#Now only the tiwi-islands manp
with nc(os.path.join(geodataf,'Topo.nc')) as topof:
    topo=np.ma.masked_equal(topof.variables['topo'][:],0)
    topolon=topof.variables['longitude'][:]
    topolat=topof.variables['latitude'][:]
with nc(CPOLF) as fnc:
    lon=fnc.variables['longitude'][:]
    lat=fnc.variables['latitude'][:]

fig = plt.figure(figsize=(20,20))
ax = fig.add_subplot(111)
ls = LightSource(azdeg=315, altdeg=45)

M = Basemap(projection='rotpole', llcrnrlat=min(Gtopolat), llcrnrlon=min(Gtopolon), urcrnrlat=max(Gtopolat), urcrnrlon=max(Gtopolon), 
            resolution='c', area_thresh=1, ax=ax, lon_0=0, o_lon_p=0., o_lat_p=90.,)
borders = (([122.963, 139.055], [-17.091, -7.443]), 
           ([127.385, 134.693], [-14.505, -10.005]), 
           ([129.557, 132.529], [-13.757, -10.757]))
text=(('UM 4km', fontsize+4), ('UM 1.33km', fontsize+4), ('UM 0.44km', fontsize+4))
M.pcolormesh(Gtopolon, Gtopolat, ls.hillshade(Gtopo[:], vert_exag=1), cmap= topocmap, ax=ax)
#M.drawmapscale(M.boundarylons[0]-0.7, M.boundarylats[0]+0.7, max(topolon), min(topolat), 10, 
#              barstyle='fancy', fontsize=19, labelstyle='simple')
for conf in range(3):
    lons,lats = borders[conf]
    txt, fts = text[conf]
    x, y = M([lons[0], lons[0], lons[1],lons[1]], [lats[0], lats[1], lats[1], lats[0]])
    xy = list(zip(x,y))
    ttt = plt.text(*M(lons[0]+0.1, lats[0]+0.05), txt, fontsize=fts)
    ttt.set_bbox(dict(facecolor='w', alpha=0.7, edgecolor='none', lw=1))
    poly = Polygon(xy, alpha=1, edgecolor='k', lw=1, facecolor='none')
    plt.gca().add_patch(poly)
#M.drawcoastlines()
M.scatter(*M(131.043101, -12.250323), marker='*',color='k', s=200)
plt.text(*M(borders[0][0][0]+0.15, borders[0][-1][1]-0.6), 'b)', fontsize=fontsize+12)
axins = zoomed_inset_axes(ax, 7, loc=2)
axins.set_xlim(min(topolon), max(topolon))
axins.set_ylim(min(topolat), max(topolat))
tmap = Basemap(projection='cyl',llcrnrlat=min(lat), llcrnrlon=min(lon), urcrnrlat=max(lat), urcrnrlon=max(lon), 
               resolution='c',  area_thresh=1, ax=axins)
tmap.pcolormesh(topolon, topolat, ls.hillshade(topo[:], vert_exag=1), cmap= topocmap)
mark_inset(ax, axins, loc1=3, loc2=4, fc="k", ec='k', lw=0.5)
plt.annotate('a)', xy=(0.01 ,0.92), xycoords='axes fraction', fontsize=fontsize+12)
lons, lats = [min(Gtopolon), max(Gtopolon)], [min(Gtopolat), max(Gtopolat)]
axins2 = inset_axes(ax, width='20%', height='20%', loc=4)
Worldmap= Basemap(projection='ortho', lat_0=Gtopolat[int(len(Gtopolat)/2)], lon_0=Gtopolon[int(len(Gtopolon)/2)],
                 ax=axins2)
Worldmap.drawcoastlines(color='k', linewidth=0.7)
Worldmap.fillcontinents(color='gray')
bx, by = Worldmap(M.boundarylons, M.boundarylats)
xy = list(zip(bx,by))
poly = Polygon(xy, alpha=1, edgecolor='none', lw=0, facecolor='darkcyan')
_ = Worldmap.ax.add_patch(poly)
_ = fig.subplots_adjust(bottom=0.35)
_ = fig.savefig('Vid/Topo.png', bbox_inches='tight', format='png', dpi=300)
fontsize = fontsize_tmp

Now calculate the map of spectral variance in the diurnal cycle

2.1 The diurnal cycle in different Monsoon stages

Burst vs break

# Define extreme Events in the CPOL era
extr1h = pd.read_hdf(extremeTS, 'P10min')
monsoon10m = pd.read_pickle(BurstF)
monsoon10m = monsoon10m.dropna()
monsoon1h = monsoon10m.groupby(pd.TimeGrouper('1h')).mean()
breaks = monsoon10m.loc[monsoon10m == 0].index
bursts = monsoon10m.loc[monsoon10m > 0].index
thresh = perc['10min'][95]
S = pd.Series(np.arange(len(times)),index=times)
print('Bursts: %0.2f %% Breaks: %0.2f %%' 
      %(len(bursts)/len(monsoon10m.index) * 100, len(breaks)/len(monsoon10m.index) * 100 ))

Bursts: 39.73 % Breaks: 60.27 %

#Plot the diurnal Cycle
%matplotlib inline
nplot=1

import matplotlib.gridspec as gridspec
gs = gridspec.GridSpec(5, 3)

sns.set_style("whitegrid", {'axes.grid':True, 'ticks':True})
from matplotlib.ticker import ScalarFormatter
#Plot the output map
fig = plt.figure(figsize=(15,8))
cmap = colmap2.Blues
cmap2 = colmap.GMT_drywet
cmap.set_under('w'),cmap.set_bad('w')
cmap2.set_under('w'),cmap2.set_bad('w')
ax1 = plt.subplot(gs[:2, 0])
#ax1 = fig.add_subplot(231)
ax1.set_title('Rainfall Climatology', size=fontsize)
X, Y = np.meshgrid(lon.filled(-50), lat.filled(-50))
im = m1.contourf(X, Y, map1, np.arange(0,13,2), cmap=colmap.GMT_drywet, ax=ax1)

#im = m1.pcolormesh(lon.filled(-50), lat.filled(-50), map1, vmin=0, cmap=colmap.GMT_drywet, ax=ax1)
m1.drawcoastlines(linewidth=.7)
cbar = m1.colorbar(im , location='bottom', pad="5%")
cbar.set_label('Rain-Rate [mm/d]', size=fontsize)
cbar.ax.tick_params(labelsize=fontsize-4)

sns.set_style("darkgrid", {'axes.grid':True, 'ticks':True})
ax2 = plt.subplot(gs[:2, 1])
sns.lineplot('x', "y", data=sns_data, ax=ax2, lw=1)
ax2.set(xscale="log")
ax2.set_xlim(0.25,60)
ax2.set_ylim(0,0.6)
ax2.set_xlabel('Periods ($\\tau$) [d]',size=fontsize)
ax2.set_title('Spectral variance [\%]',size=fontsize)
ax2.tick_params(labelsize=fontsize-4)
ax2.set_ylabel('')
ax2.xaxis.set_major_formatter(ScalarFormatter())
sns.set_style("whitegrid", {'axes.grid':True, 'ticks':True})

ax3 = plt.subplot(gs[:2, 2])
im2 = m1.contourf(X, Y, map2, np.arange(0, 9)/10. ,cmap='Blues', ax=ax3)
m1.drawcoastlines(linewidth=.7)
ax3.set_title('Diurnal Cycle', size=fontsize)
cbar = m1.colorbar(im2,location='bottom',pad="5%")
cbar.set_label('Spectral variance [\%]', size=fontsize)
cbar.ax.tick_params(labelsize=fontsize-4)
# %%
#fig.subplots_adjust(right=0.98, bottom=0.03, top=0.99,left=0.01, hspace=0, wspace=0.17)
#'''
sns.set_style("darkgrid", {'axes.grid':True, 'ticks':True})
#fig, ax = plt.subplots(1, figsize=(15,10))
ax4 = plt.subplot(gs[3:, :])
sns.lineplot(x='Local Time', y='Rain-Rate', data=diurnacl_cycle, hue='Monsoon-Stage', lw=3, ax=ax4)
ax4.set_ylabel('Rain-Rate mm/d', fontsize=fontsize)
ax4.set_xlabel('Local Time [h]', fontsize=fontsize)
#ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%_H'))
ax4.legend(loc=0, fontsize=fontsize-4)
ax4.tick_params(labelsize=fontsize)
fig.subplots_adjust(bottom=0.01, top=0.98, left=0.01, right=0.98, hspace=0.1, wspace=0.2)
#fig.subplots_adjust(bottom=0.145, top=0.99, right=0.99, left=0.1)
nplot += 1; _ = fig.savefig('Vid/Fig_%03i.png'%nplot, bbox_inches='tight', format='png', dpi=300)

2.2 Local Extemes and the stage of the Monsoon

The above plots show that diurnal cycle is the most important mode of variability in the considered area. The island thunderstrom hector is the most extreme convective system that occurs at the time-scale at perdiods of 24 h. The Hector 'signal' is most pronounced during the Monsoon break period. Before studying Hector we split the CPOL area into break and burst period and investigate the occurrence of extremes during Monsoon break and burst.

First we define extreme events by applying the 99th percentile threshold of the total 10min dataset. Then we define the break and burst events by the definition of Narsey et al. and calculate the number of burst and the number of break events during the CPOL era

# Now Plot the maps
mpld3.disable_notebook()
sns.set_style("whitegrid", {'axes.grid':True, 'ticks':True})
maxval =[(0,15), (0,15), (-12,12), (0, 150), (0, 150), (-200, 200)]
names=['Breaks', 'Bursts', 'Breaks - Bursts',
      'Extremes During Breaks', 'Extremes During Bursts', 'Breaks - Bursts']
cmap_list = [(colmap.GMT_drywet,'k'), (colmap.GMT_drywet,'k'), (colmap.GMT_polar_r,'k'),
       (mpl.cm.magma,'w'), (mpl.cm.magma,'w'), (colmap.GMT_polar_r,'k')]

C = [Mean[1]*1.5, Mean[2], Mean[1]*1.5 - Mean[2], D[1]*1.2, D[2], D[1]*1.2 - D[2]]
fig = plt.figure(figsize=(15,10))
mask = ((D[0] * 0 )+1)
im=[]
for i, idx in enumerate(C):
    ax = fig.add_subplot(2,3,i+1)
    data = np.zeros([len(lat),len(lon)])
    colm, coast = cmap_list[i]
    colm.set_bad(dict(k='w', w='k')[coast])
    colm.set_under(dict(k='w', w='k')[coast])
    try:
        im.append(m.pcolormesh(lon, lat, mask*C[i].filled(0),vmin=maxval[i][0],vmax=maxval[i][1],cmap=colm))
    except AttributeError:
        im.append(m.pcolormesh(lon, lat, mask*C[i],vmin=maxval[i][0],vmax=maxval[i][1], cmap=colm))
    m.drawcoastlines(color=coast)
    ax.set_title('%s'%names[i],size=fontsize)
fig.subplots_adjust(right=0.98, bottom=0.01, top=0.98,left=0.01, hspace=0, wspace=0)
cbar_ax1 = fig.add_axes([0.05, 0.57,0.58, 0.015])
cbar1 = fig.colorbar(im[0], cax=cbar_ax1, orientation='horizontal')
cbar1.ax.tick_params(labelsize=fontsize-4)
cbar1.set_label(u'Rain-Rate [mm/d]',size=fontsize)
cbar_ax2 = fig.add_axes([0.685, 0.57,0.27, 0.015])
cbar2 = fig.colorbar(im[2], cax=cbar_ax2, orientation='horizontal')
cbar2.ax.tick_params(labelsize=fontsize-4)
cbar2.set_label(u'Difference',size=fontsize)
cbar_ax3 = fig.add_axes([0.05, 0.09,0.58, 0.015])
cbar3 = fig.colorbar(im[3], cax=cbar_ax3, orientation='horizontal')
cbar3.ax.tick_params(labelsize=fontsize-4)
cbar3.set_label(u'N. of  events above 99th%',size=fontsize)
cbar_ax4 = fig.add_axes([0.685, 0.09,0.27, 0.015])
cbar4 = fig.colorbar(im[-1], cax=cbar_ax4, orientation='horizontal')
cbar4.ax.tick_params(labelsize=fontsize-4)
cbar4.set_ticks([-200, -100, 0, 100, 200])
cbar4.set_label(u'Difference',size=fontsize)
nplot+=1 ; _ = fig.savefig('Vid/Fig_%03i.png'%nplot, bbox_inches='tight', format='png', dpi=300)

2.1 A Week of Hector

Now we compare the ensemble timeseries of the simulated rainfall for the 1.33km and 0.44km Simulations with the CPOL observations.

time, obs, rain = [], [], []
for e in range(len(ensembles)):
    r1 = list(UM133ens.variables['lsrain'][e].mean(axis=(1,2,3)).values)
    r2 = list(UM044ens.variables['lsrain'][e].mean(axis=(1,2,3)).values)
    obs+=len(r1)*['UM 1.33km'] + len(r2)*['UM 0.44km']
    t1 = pd.DatetimeIndex(UM133ens.coords['t'].values).tz_localize(utc).tz_convert(timezone).tz_localize(None).to_pydatetime()
    rain += r1 + r2
    t2 = pd.DatetimeIndex(UM044ens.coords['t'].values).tz_localize(utc).tz_convert(timezone).tz_localize(None).to_pydatetime()
    time += list(t1) + list(t2)

r3 = list(OBS.dataset[1].variables['lsrain'].mean(axis=(1,2)).values)
t3 = list(pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).tz_localize(None).to_pydatetime())
o3 = len(r3)*['CPOL']
rain += r3
obs += o3
time += t3
plot_data = pd.DataFrame({'rain-rate':rain, 'time':time, 'Type':obs})

#Plot area avg TS
import matplotlib.dates as mdates
mpld3.disable_notebook()
sns.set_style("darkgrid", {'axes.grid':True, 'ticks':True})
myFmt = mdates.DateFormatter('%d/%m')
#
#Get the minimum time period that is covered by all simulations and also observations
fig = plt.figure(figsize=(15,7))
ax = fig.add_subplot(111)
sns.lineplot(x='time', y='rain-rate', hue='Type', ax=ax, data=plot_data, lw=1.3)

plt.subplots_adjust(left=0.1, bottom=0.1, right=0.98, top=0.9)
ax.legend(loc=0, fontsize=fontsize-4)
ax.tick_params(labelsize=fontsize-4)
#ax.set_xlim(plot_data.time.iloc[0], plot_data.time.iloc[-1])
ax.set_ylim(-0.001,0.5)
ax.set_ylabel('Rain-Rate [mm/h]', fontsize=fontsize)
_ = ax.xaxis.set_major_formatter(myFmt)
_ = ax.set_xlabel('')
#_ = ax.set_title('Comparison of Area-Avg Rainfall', fontsize=fontsize)
fig.subplots_adjust(bottom=0.145, top=0.9, right=0.99, left=0.1)
nplot+=1; _ = fig.savefig('Vid/Fig_%03i.png'%nplot, bbox_inches='tight', format='png', dpi=300)

2.1.1 Maps of Rainfall

Now lets look at some animations of UM rainfall occuring over the Tiwi-islands and compare it to the cpol observations

os.chdir('/home/unimelb.edu.au/mbergemann/Work/Programming/Extremes/python/')

#Plot Animatoin
import io
import base64
from IPython.display import HTML
outvid=os.path.join('Vid','WeekOfHector-Ens-3.mp4')
video = io.open(outvid, 'r+b').read()
encoded = base64.b64encode(video)
HTML(data='''<video alt="test" width="950" height="500" loop="true" controls>
                <source src="data:video/mp4;base64,{0}" type="video/mp4" />
             </video>'''.format(encoded.decode('ascii')))

Maps of Rainfall (0.44km)

#Create the avg. maps of rainfall 0.44km
#m1, m2 = None, None
mpld3.disable_notebook()
sns.set_style("whitegrid", {'axes.grid':True, 'ticks':True})
rom = {1:'i', 2:'ii', 3:'iii', 4:'iv', 5:'v', 6:'vi', 7:'vii', 8:'viii'}
fig = plt.figure(figsize=(15,12))

UM = UM044ens.groupby('t.day').sum(dim='t').mean(dim='day')
obs = OBS.dataset[1].groupby('t.day').sum(dim='t').mean(dim='day')
colm = colmap2.Blues

colm.set_under('w', alpha=0)
colm.set_bad('w', alpha=0)
tsteps = pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).to_pydatetime()
obs_data = np.ma.masked_less(obs.variables['lsrain'][:].values, -0.01) / 6
um_data = np.ma.masked_less(UM.variables['lsrain'][:].values, -0.01) / 6
lat1, lon1 = obs.variables['latitude'][:,0], obs.variables['longitude'][0,:]
lat2, lon2 = UM.coords['lat'][:], UM.coords['lon'][:]
ax = [fig.add_subplot(3,3,1)]
if m2 is None:
    m2 = Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2),
             resolution='f', area_thresh=1)
if m1 is None:
    m1 = Basemap(llcrnrlat=min(lat1), llcrnrlon=min(lon1), urcrnrlat=max(lat1), urcrnrlon=max(lon1),
             resolution='f', area_thresh=1)
#m[0].pcolormesh(topolon, topolat, ls.hillshade(topo[:], vert_exag=1), cmap='gray', ax=ax[0])
im = [m2.pcolormesh(lon1, lat1, obs_data,vmin=0.0,vmax=5,cmap=colm)]
m2.drawcoastlines(linewidth=0.8)
cbar_ax = fig.add_axes([0.12, 0.12, 0.74, 0.02])
cbar=fig.colorbar(im[-1], cax=cbar_ax, orientation='horizontal')
cbar.ax.tick_params(labelsize=24)
cbar.set_label('Rain-Rate [mm/d]',size=24)
tit = ['CPOL']
ax[-1].set_title('CPOL', fontsize=fontsize)
for i in range(1,9):
    if m2 is None:
        m2 = Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2),
                         resolution='f', area_thresh=1)
    ax.append(fig.add_subplot(3,3,i+1))
    m2.drawcoastlines(linewidth=0.8)
    im.append(m2.pcolormesh(lon2, lat2, um_data[i-1][0],vmin=0.0,vmax=5,cmap=colm))
    tit.append('Ens. %s'%rom[i].upper())
    ax[-1].set_title(tit[-1],fontsize=fontsize)
fig.subplots_adjust(right=0.99, bottom=0.15, top=0.95,left=0.01, hspace=0.2, wspace=0.0)
nplot+=1 ; _ = fig.savefig('Vid/Fig_%03i.png'%nplot, bbox_inches='tight', format='png', dpi=300)

Maps of Rainfall (1.33km)

#Create the avg maps of rainfall for 1.33km
sns.set_style("whitegrid", {'axes.grid':True, 'ticks':True})
fig = plt.figure(figsize=(15,12))
UM = UM133ens.groupby('t.day').sum(dim='t').mean(dim='day')
obs = OBS.dataset[1].groupby('t.day').sum(dim='t').mean(dim='day')
colm = colmap2.Blues
colm.set_under('w', alpha=0)
colm.set_bad('w', alpha=0)
tsteps = pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).to_pydatetime()
obs_data = np.ma.masked_less(obs.variables['lsrain'][:].values, -0.01) / 6
um_data = np.ma.masked_less(UM.variables['lsrain'][:].values, -0.01) / 6
lat1, lon1 = obs.variables['latitude'][:,0], obs.variables['longitude'][0,:]
lat2, lon2 = UM.coords['lat'][:], UM.coords['lon'][:]
ax = [fig.add_subplot(3,3,1)]
if m1 is None:
    m1 = Basemap(llcrnrlat=min(lat1), llcrnrlon=min(lon1), urcrnrlat=max(lat1), urcrnrlon=max(lon1), resolution='f',
                 area_thresh=1)
#m[0].pcolormesh(topolon, topolat, ls.hillshade(topo[:], vert_exag=1), cmap='gray', ax=ax[0])
im = [m2.pcolormesh(lon1, lat1, obs_data,vmin=0.0,vmax=5,cmap=colm)]
m2.drawcoastlines(linewidth=0.8)
cbar_ax = fig.add_axes([0.12, 0.12, 0.74, 0.02])
cbar=fig.colorbar(im[-1], cax=cbar_ax, orientation='horizontal')
cbar.ax.tick_params(labelsize=24)
cbar.set_label('Avg. Rain-Rate [mm/d]',size=24)
tit = ['CPOL']
ax[-1].set_title('CPOL', fontsize=fontsize)
for i in range(1,9):
    if m2 is None:
        m2 = Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2), resolution='f',
                     area_thresh=1)
    ax.append(fig.add_subplot(3,3,i+1))
   
    m2.drawcoastlines(linewidth=0.8)
            #m[-1].shadedrelief()
    im.append(m2.pcolormesh(lon2, lat2, um_data[i-1][0],vmin=0.0,vmax=5,cmap=colm))
    tit.append('Ens. %s'%rom[i].upper())
    ax[-1].set_title(tit[-1],fontsize=fontsize)
fig.subplots_adjust(right=0.99, bottom=0.15, top=0.95,left=0.01, hspace=0.2, wspace=0.0)
nplot+=1 ; _ = fig.savefig('Vid/Fig_%03i.png'%nplot, bbox_inches='tight', format='png', dpi=300)

Ensemble mean and std

#Plet the avg maps for comparison
fig = plt.figure(figsize=(18,12))
#sns.set_style("darkgrid", {'axes.grid':True, 'ticks':True})
UM_2 = UM133ens.groupby('t.day').sum(dim='t').mean(dim='day')
UM_1 = UM044ens.groupby('t.day').sum(dim='t').mean(dim='day')
obs = OBS.dataset[1].groupby('t.day').sum(dim='t').mean(dim='day')
colm = colmap2.Blues
colm.set_under('w', alpha=0)
tsteps = pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).to_pydatetime()
obs_data = np.ma.masked_less(obs.variables['lsrain'][:].values, -0.01)
data = [obs_data[1:-1]/6,
        np.ma.masked_less(UM_1.variables['lsrain'][:].values, -0.01).mean(axis=0)[0]/6,
        np.ma.masked_less(UM_2.variables['lsrain'][:].values, -0.01).mean(axis=0)[0]/6,
        None,
        np.ma.masked_less(UM_1.variables['lsrain'][:].values/6, -0.01).std(axis=0)[0],
        np.ma.masked_less(UM_2.variables['lsrain'][:].values/6, -0.01).std(axis=0)[0]
       ]

data.append(None)

data.append(data[0] - data[1])
data.append(data[0] - data[2])
UM_1min = np.ma.masked_less(UM_1.variables['lsrain'][:].values, -0.01).min(axis=0)[0]
UM_2min = np.ma.masked_less(UM_2.variables['lsrain'][:].values, -0.01).min(axis=0)[0]

UM_1max = np.ma.masked_less(UM_1.variables['lsrain'][:].values, -0.01).max(axis=0)[0]
UM_2max = np.ma.masked_less(UM_2.variables['lsrain'][:].values, -0.01).max(axis=0)[0]

#data.append(None)
#data.append(UM_1max - UM_1min)
#data.append(UM_1max - UM_1min)


lat1, lon1 = obs.variables['latitude'][:,0], obs.variables['longitude'][0,:]
lat2, lon2 = UM_1.coords['lat'][:], UM_1.coords['lon'][:]
ax = [fig.add_subplot(3,3,1)]
#m[0].pcolormesh(topolon, topolat, ls.hillshade(topo[:], vert_exag=1), cmap='gray', ax=ax[0])
im = [m2.pcolormesh(lon1, lat1, obs_data/6,vmin=0.0,vmax=25/6,cmap=colm)]
m2.drawcoastlines()
#cbar_ax = fig.add_axes([0.14, 0.33, 0.74, 0.01])
#cbar=fig.colorbar(im[-1], cax=cbar_ax, orientation='horizontal')
#cbar.ax.tick_params(labelsize=24)
#cbar.set_label('Avg. Rain-rate [mm/d]',size=24)
tit = ['CPOL avg.', 'UM 0.44km ens avg.', 'UM 1.33km ens avg.',
       None, 'UM 0.44 ens std.', 'UM 1.33km ens std.',
       None, 'CPOL - UM 0.44 ens avg.' ,'CPOL - UM 1.33km ens avg.',
       None, 'UM 0.44 ens range', 'UM 1.33km ens range']
colors = {1:colm, 2:colm, 4:colm, 5:colm,
          7:colmap.GMT_polar_r, 8:colmap.GMT_polar_r,
         10:colm, 11:colm}
Range = {1:(0.0, 25/6), 2:(0.0, 25/6), 4:(0,8/6), 5:(0,8/6), 7:(-20/6,20/6), 8:(-20/6,20/6), 10:(0,30/6), 11:(0,30/6)}
height = {2:0.682, 5:0.37, 8:0.044, 11:0.05}
ax[-1].set_title('CPOL avg.', fontsize=fontsize-2)
cbar_ax = [fig.add_axes([0.91, height[2], 0.01, 0.23])]
cbar = [fig.colorbar(im[-1], cax=cbar_ax[-1], orientation='vertical')]
cbar[-1].ax.tick_params(labelsize=fontsize-4)
cbar[-1].set_label('Rain-Rate [mm/d]',size=fontsize-4)
for i in range(1,len(data)):
    if data[i] is not None:
        ax.append(fig.add_subplot(3,3,i+1))
            #m[-1].shadedrelief()
        im.append(m2.pcolormesh(lon2, lat2, data[i],vmin=Range[i][0],vmax=Range[i][1],cmap=colors[i]))
        m2.drawcoastlines()
        ax[-1].set_title(tit[i],fontsize=fontsize-2)
        if i in (5, 8, 11):
            cbar_ax.append(fig.add_axes([0.91, height[i], 0.01, 0.23]))
            cbar.append(fig.colorbar(im[-1], cax=cbar_ax[-1], orientation='vertical'))
            cbar[-1].ax.tick_params(labelsize=fontsize-4)
            cbar[-1].set_label('Rain-rate [mm/d]',size=fontsize-2)
fig.subplots_adjust(right=0.9, bottom=0.01, top=0.95,left=0.01, hspace=0.0, wspace=0)
nplot+=1 ; _ = fig.savefig('Vid/Fig_%i03.png'%nplot, bbox_inches='tight', format='png', dpi=300)

2.1.2 The Simulated Diurnal Cycle

data, time, model = [], [], []
sns.set_style("darkgrid", {'axes.grid':True, 'ticks':True})
for i in range(len(ensembles)):
    d1 = UM133.dataset[i].groupby('t.hour').mean(dim='t').mean(dim=('surface','lat','lon'))
    d2 = UM044.dataset[i].groupby('t.hour').mean(dim='t').mean(dim=('surface','lat','lon'))
    t1 = pd.Series(d1.variables['lsrain'].values, index=(d1.variables['hour'].values+9) % 24 ).sort_index(inplace=False)
    t2 = pd.Series(d2.variables['lsrain'].values, index=(d2.variables['hour'].values+9) % 24 ).sort_index(inplace=False)
    time += list(t1.index) + list(t2.index)
    data += list(t1.values) + list(t2.values)
    model += len(t1)*['UM 1.33km'] + len(t2)*['UM 0.44km']
cpol_1 = OBS.dataset[1].groupby('t.hour').mean(dim='t').mean(dim=('y','x'))
obs_S = pd.Series(cpol_1.variables['lsrain'].values, index=(cpol_1.variables['hour'].values+9) % 24 ).sort_index(inplace=False)
data += list(obs_S.values)
time += list(obs_S.index)
model += len(obs_S)*['Cpol']
plot_data = pd.DataFrame(dict(data=data, time=time, Type=model))
#nplot+=1 ; _ = fig.savefig('Vid/Fig_%03i.png'%nplot, bbox_inches='tight', format='png', dpi=300)
#mpld3.enable_notebook()

Location and timing of daily rainfall peaks

#Plot the maps
plt.close()
mpld3.disable_notebook()
sns.set_style("whitegrid", {'axes.grid':True, 'ticks':False})
um_1 = UM044ens.groupby('t.hour').mean(dim='t').mean(dim=('surface'))
um_2 = UM133ens.groupby('t.hour').mean(dim='t').mean(dim=('surface'))
obs  = OBS.dataset[1].groupby('t.hour').mean(dim='t')

time = (um_1.hour + 9) % 24

data = [obs, um_2, um_1]

fig = plt.figure(figsize=(15,7))
lat1, lon1 = OBS.dataset[1].variables['latitude'].values[:,0], OBS.dataset[1].variables['longitude'].values[0,:]
lat2, lon2 = um_1.coords['lat'][:], um_1.coords['lon'][:]
fig.subplots_adjust(right=0.94, bottom=0.025, top=0.98,left=0.01, hspace=0, wspace=0)
first = True
val_max = 2
colm = colmap2.Blues
colm.set_under('w', alpha=0)
colm.set_bad('w', alpha=0)

#ax = [fig.add_subplot(2,3,1)]
mm, im = [], []
col = 2
for i in range(24):
    #break
    fname_1 = os.path.join(outdir, 'DiurnalCycle_%02i.png'%i)
    if first:
        for col in range(1, 7):
            if col != 4:
                ax1 = fig.add_subplot(2,3,col)
                m2.drawcoastlines(linewidth=0.8)
            if col == 1:
                ax1.set_title('CPOL',fontsize=fontsize)
                _ = im.append(m2.pcolormesh(lon1, lat1, obs['lsrain'].values[i],  shading='gouraud', vmin=0.1, vmax=val_max, cmap=colm))
            elif col == 2:
                ax1.set_title('UM 1.33km mean', fontsize=fontsize)
                _ = im.append(m2.pcolormesh(lon2, lat2, um_2['lsrain'].mean(dim='ens').values[i],  shading='gouraud', vmin=0.1, vmax=val_max, cmap=colm))
            elif col == 3:
                ax1.set_title('UM 0.44km mean', fontsize=fontsize)
                _ = im.append(m2.pcolormesh(lon2, lat2, um_1['lsrain'].mean(dim='ens').values[i],  shading='gouraud', vmin=0.1, vmax=val_max, cmap=colm))
            elif col == 5:
                ax1.set_title('UM 1.33km std', fontsize=fontsize)
                _ = im.append(m2.pcolormesh(lon2, lat2, um_2['lsrain'].std(dim='ens').values[i],  shading='gouraud', vmin=0.1, vmax=val_max, cmap=colm))
            elif col == 6:
                ax1.set_title('UM 0.44km std', fontsize=fontsize)
                _ = im.append(m2.pcolormesh(lon2, lat2, um_2['lsrain'].std(dim='ens').values[i],  shading='gouraud', vmin=0.1, vmax=val_max, cmap=colm))      
         
                cbar_ax = fig.add_axes([0.14, 0.0, 0.74, 0.02])
                cbar=fig.colorbar(im[-1], cax=cbar_ax, orientation='horizontal')
                _ = cbar.ax.tick_params(labelsize=fontsize-4)
        first = False
    else:
        ###MEAN
        im[0].set_array(obs['lsrain'].values[i].ravel())
        im[1].set_array(um_2['lsrain'].mean(dim='ens').values[i].ravel())
        im[2].set_array(um_1['lsrain'].mean(dim='ens').values[i].ravel())
        ###STD
        im[-2].set_array(um_2['lsrain'].std(dim='ens').values[i].ravel())
        im[-1].set_array(um_1['lsrain'].std(dim='ens').values[i].ravel())
    cbar.set_label('Rain-Rate at %2i LT [mm/h]'%((i+9.5)%24),size=fontsize)
    fig.savefig(fname_1, bbox_inches='tight', format='png', dpi=72)
    #break
dest_dir = os.path.abspath('Vid')
make_mp4_from_frames(outdir, dest_dir, 'WeekOfHector-Diurnal-2.mp4', 4, glob='DiurnalCycle_??')
fig.clf()
plt.close()
day = (time[1:5] - 9) % 24
night = (time[5:5+4] - 9) % 24
obs_diff = obs['lsrain'][day].sum(dim='hour') - obs['lsrain'][night].sum('hour')
um_1_diff = um_1['lsrain'][day].sum(dim='hour') - um_1['lsrain'][night].sum('hour')
um_2_diff = um_2['lsrain'][day].sum(dim='hour') - um_2['lsrain'][night].sum('hour')

mask = np.ma.masked_invalid(obs['lsrain'].values) * 0 + 1
obs_max = (((np.argmax(obs['lsrain'].values * mask, axis=0)*mask[0] + 9.5 ) % 24)/2).round(0) * 2
um_1_max = ((um_1['lsrain'].argmax(dim='hour').mean(dim='ens').round(0) + 9.5) % 24 / 2).round(0) * 2 
um_2_max = ((um_2['lsrain'].argmax(dim='hour').mean(dim='ens').round(0) + 9.5) % 24 / 2).round(0) * 2
colm = CubeYF_20.get_mpl_colormap(N=24, gamma=2.0)
#colm = qualitative.Paired[12].get_mpl_colormap(N=24, gamma=2.0)
colm.set_bad('w')
#m, ax = [], []
ax = []
data = [obs_max, um_1_max, um_2_max]
names = ['CPOL', 'UM 1.33km', 'UM 0.44km']
gs = gridspec.GridSpec(30, 3)

fig = plt.figure(figsize=(15,8))
fig.subplots_adjust(right=0.98, bottom=0.03, top=0.98,left=0.01, hspace=0, wspace=0)
cbar_ax = plt.subplot(gs[13, :])

for i in range(3):
    ax.append(plt.subplot(gs[:13, i]))
    if m2 is None:
        m2 = Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2), resolution='f',
                     area_thresh=1)
    
    m2.drawcoastlines()
    try:
        im = m2.pcolormesh(lon1, lat1, data[i], vmin=0, vmax=24, cmap=colm)
    except TypeError:
        im = m2.pcolormesh(lon2, lat2, data[i], vmin=0, vmax=24, cmap=colm)
    ax[-1].set_title(names[i], fontsize=fontsize)
cbar=fig.colorbar(im, cax=cbar_ax, orientation='horizontal',ticks=list(range(1,24,2)))
cbar.ax.tick_params(labelsize=fontsize-4)
#cbar.set_label('Local Time [h]', size=fontsize)
mpld3.disable_notebook()
# Plot avg diurnal cycle
#Get the minimum time period that is covered by all simulations and also observations
#fig = plt.figure(figsize=(15,8), dpi=72)
#ax = fig.add_subplot(111)
sns.set_style("darkgrid", {'axes.grid':True, 'ticks':True})
ax = plt.subplot(gs[18:, :])

####################
sns.lineplot(x=plot_data.time, y=plot_data.data, hue=plot_data.Type, ax=ax, lw=3, data=plot_data)
#ax.set_ylim(0.0,0.42)
ax.set_xlim(0,23)
ax.set_ylabel('Rain-Rate [mm/h]', fontsize=fontsize)
ax.set_xlabel('Local Time [h]', fontsize=fontsize)
ax.legend(loc=0, fontsize=fontsize-4)
ax.tick_params(labelsize=fontsize-4)
fig.subplots_adjust(bottom=0.145, top=0.99, right=0.99, left=0.1)
#ax.grid(color='r', linestyle='-', linewidth=0)
mpld3.disable_notebook()
nplot = 12; _ = fig.savefig('Vid/Fig_%03i.png'%nplot, bbox_inches='tight', format='png', dpi=300)

3.1 Pdf's of Rainfall

#Plot the pdfs
sns.set_style('darkgrid')
from mpl_toolkits.axes_grid.inset_locator import inset_axes
dist=lognorm([LogNorm.CPOL['shape']*3.4], loc = np.log(LogNorm.CPOL['scale']*1.9))#, shape  # mu, sigma
med = np.median(um044vec.loc[um044vec>0].values)
X = np.linspace(0,65,200)
#b = 2.5
nbins = 100
def fitpdf(x, b):
    a=med*100
    return ((b/a) * (x/a)**(b-1)) / (1+(x/a)**b)**2
#popt, pcov = curve_fit(powerlaw, um044cnt['Rain-rate'].values[1:],
#                       um044cnt['count'].values[1:])
popt, pcov = PowerFit.CPOL['popt']/4, PowerFit.CPOL['pcov']
fig = plt.figure(figsize=(15,10))
ax = fig.add_subplot(111)
ax2 = ax.twinx()
#Define histogram data to plot
histdata = (um133vec.loc[um133vec>0], um044vec.loc[um044vec>0], obs_vec.loc[obs_vec>0])
#facecolors = ['fuchsia', 'firebrick', 'lightseagreen']
#n1_1, bins1_1, patches1_1 = ax2.hist(, bins=nbins, normed=1, facecolor='firebrick', alpha=0.25)
#n2, bins2, patches2 = ax2.hist(, bins=nbins, normed=1, facecolor='lightseagreen', alpha=0.25)
ax2.set_ylim(0,0.9)
ax2.set_yticks([])
#ax.plot(np.linspace(0.01, 100,200), fitpdf(np.linspace(0.01,100,200),popt[0]), lw=4, color='navy',
#        label='Log-logistic $\\beta = %02.2e$'%popt[0])
#ax.plot(cnts.UM044['Rain-rate'].values, cnts.UM044['counts'].values,'-', label='UM 0.44 km ens',lw=4, color='fuchsia')
#ax.plot(cnts.UM133['Rain-rate'].values, cnts.UM133['counts'].values,'--', label='UM 1.33 km ens',lw=4, color='fuchsia')
#ax.plot(cnts.CPOL['Rain-rate'].values, cnts.CPOL['counts'].values,'-', label='CPOL',lw=4, color='lightseagreen')
#ax.scatter(um044cnt['Rain-rate'].values, um044cnt['count'].values, c='red', s=125, alpha=0.24)
#ax.scatter(obscnt['Rain-rate'].values, obscnt['count'].values, c='purple', s=125, alpha=0.24)
#ax.plot(np.linspace(0, 50,200), um044pdf, lw=4, label='UM 0.44km ens', color='lightseagreen')
#ax.plot(np.linspace(0, 50,200), obspdf, lw=4, label='CPOL', color='fuchsia')
#ax.plot(np.linspace(0.0,50,200), dist.pdf(np.linspace(0.0, 50,200)), color='firebrick', alpha=0.8,
#        lw=4, label='Log-norm ($\\mu=%02.2f, \\sigma=%02.2f$)'%(scale, shape))
n1, bins1, patches1 = ax.hist(histdata, bins=nbins, normed=1, alpha=0.55)
#fitted_pl.plot_pdf(ax=ax,color='k',lw=4)
#ax.plot(np.linspace(0.1, 100), um044dens_aavg(np.linspace(0.1,100)), lw=4, label='UM 0.44km area-avg')
#ax.plot(np.linspace(0.1, 100), um133dens_aavg(np.linspace(0.1,100)), lw=4, label='UM 1.33km area-avg')


ax.tick_params(labelsize=fontsize-4)

#ax.set_ylim(0.,.25)
ax.set_xlim(0.,20)
ax.set_ylabel('Density', fontsize=fontsize)
ax.set_xlabel('Rain-Rate [mm/h]', fontsize=fontsize)
axin = inset_axes(ax, width = "80%", height = "80%", loc=1)

axin.plot(cnts.UM133['Rain-rate'].values, cnts.UM133['counts'].values,'-', label='UM 1.33 km ens',lw=2)
axin.plot(cnts.UM044['Rain-rate'].values, cnts.UM044['counts'].values,'-', label='UM 0.44 km ens',lw=2)
axin.plot(cnts.CPOL['Rain-rate'].values, cnts.CPOL['counts'].values,'-', label='CPOL',lw=2)
#axin.plot(np.linspace(0.1,80,200), dist.pdf(np.linspace(0.1, 80,200)), lw=1, alpha=0.8,
#          label='Log-norm ($\\mu=%02.2f, \\sigma=%02.2f$)'%(scale, shape))
#axin.set_xlabel('$\\log$(Rain-rate)', fontsize=fontsize-2)
#axin.plot(X, um044pdf,'-', label='UM 0.44 km ens',lw=4, color='fuchsia')
#axin.plot(X, obspdf,'-', label='CPOL',lw=4, color='lightseagreen')

axin.set_yscale('log')
axin.set_xscale('log')
#axin.set_xticks([])
#axin.set_yticks([])
axin.tick_params(labelsize=fontsize-4)
axin.set_xlim(0.1, 100)
axin.set_ylim(0.000023, 1)
ann = axin.annotate("(b)", (0.95,0.93), xycoords="axes fraction", fontsize=fontsize)
axin.legend(loc=3, fontsize=fontsize-4)
#n.set_ylabel('$\\log$(Density)', fontsize=fontsize-2)

fig.subplots_adjust(bottom=0.13, right=0.98, top=0.97)
nplot+=1 ; _ = fig.savefig('Vid/Fig_%03i.png'%nplot, bbox_inches='tight', format='png', dpi=300)
#mpld3.enable_notebook()

3.2 Percentiles

Now get the percentiles for the ensemble simulation and compare it against the observations

#Plot the percentages
fig = plt.figure(figsize=(15,10))
ax = fig.add_subplot(111)
plt.subplots_adjust(bottom=0.05, right=0.9, top=0.4)

ax.plot(PERC.index, PERC['UM 1.33km'].values, label='UM 1.33km',lw=3)
ax.plot(PERC.index, PERC['UM 0.44km'].values, label='UM 0.44km', lw=3)
ax.plot(PERC.index, PERC['Obs'].values, label='CPOL',lw=2)
ax.set_xlim(1,100)
#ax.set_ylim(10,100)
#ax.set_yscale('log')
#ax.set_xscale('log')
ax.set_xlabel('Percentile []', fontsize=fontsize)
ax.set_ylabel('Rain-rate [mm/h]', fontsize=fontsize)

ax.tick_params(labelsize=fontsize-4)
axin = inset_axes(ax, width = "80%", height = "60%", loc=9)

axin.plot(PERC.index[1:], PERC['UM 1.33km'].values[1:], label='UM 1.33km',lw=3)
axin.plot(PERC.index[1:], PERC['UM 0.44km'].values[1:], label='UM 0.44km', lw=3)
axin.plot(PERC.index[1:], PERC['Obs'].values[1:], label='CPOL',lw=2)

#axin.set_xticks([])
#axin.set_yticks([])
axin.tick_params(labelsize=fontsize-4)
axin.set_xlim(50, 100)
axin.set_ylim(.25,100)
#axin.set_xticks([20, 200, 500])
axin.legend(loc=0, fontsize=fontsize-4)
#axin.get_xaxis().set_major_formatter(mpl.ticker.ScalarFormatter(useMathText=False))
#axin.get_xaxis().get_major_formatter().labelOnlyBase = True
#axin.set_xticklabels([], rotation=0)
ann = axin.annotate("(b)", (0.95,0.93), xycoords="axes fraction", fontsize=fontsize)
#axin.set_xticks(range(50,100,10), minor=False)
axin.set_yscale('log')
axin.set_xscale('log')
fig.subplots_adjust(bottom=0.13, right=0.98, top=0.97)
nplot+=1 ; _ = fig.savefig('Vid/Fig_%03i.png'%nplot, bbox_inches='tight', format='png', dpi=300)

nplot

15

mpl.ticker.ScalarFormatter?

2.1 Distribution of Percentiles during break and bursts

Now we compare the distributions of percentiles during break and burst periods. We also add

Created with Jupyter, delivered by Fastly, rendered by Rackspace.